An orbital-free self-consistent field approach for molecular clusters and liquids 
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We present an "orbital" free density functional theory for computing the quantum ground state of 
atomic clusters and liquids. Our approach combines the Bohm hydrodynamical description of quan- 
tum mechanics with an information theoretical approach to determine an optimal quantum density 
function in terms of density approximates to a statistical sample. The ideas of Bayesian statistical 
analysis and an expectation-maximization procedure are combined to develop approximations to 
the quantum density and thus find the approximate quantum force. The quantum force is then 
combined with a Lennard- Jones force to simulate clusters of Argon atoms and to obtain the ground 
state configurations and energies. As demonstration of the utility and flexibility of the approach, 
we compute the lowest energy structures for small rare-glass clusters. Extensions to many atom 
systems is straightforward. 
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I. INTRODUCTION 

It has been long recognized that computational effort 
of grid-based quantum mechanical methods for nuclear 
dynamical problems grows exponentially with the num- 
ber of degrees of freedom. This limits the size of systems 
that can be handled in an exact manner to those with 
4 or less atoms. This is perhaps best illustrated in the 
field of reactive scattering experiments which have been 
limited to systems with 6D 1,2,3,4 , but it is also clearly 
seen in other areas such as photodissociation processes 
In light of this, considerable progress has been made 
in developing rigorous approaches for contracting the 
basis size required to perform such calculations. One 
such approach that has seen considerable success is the 
multi-configurational time-dependent Hartree approach 
(MCTDH) developed by Meyer and co-workers^ that 
overcomes this limitation in a numerically exact way by 
expanding the time-dependent wave function in terms of 
a number of time-dependent configurations. 



*(*) = 5>.K*)M*). 



in which the single particle (or quasi-particle) basis func- 
tions and the expansion coefficients are are coupled 
by the MCTDH equations of motion. 

For condensed phase systems and liquids path inte- 
gral Monte Carlo (PIMC) and centroid based molecu- 
lar dyamics remain the method of choice. They have 
been extremely successful in calculating a wide variety of 
different thermodynamic properties of heavily quantum 



systems. 
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Despite the success of PIMC approaches 



of late there are some inherent difficulties it faces. For 
instance, at low temperatures the amount of parameters 
that must be included can become prohibitive and lead 
to slow convergence. 

We present an approach which can model low tem- 
perature Lennard-Jones clusters with ease. The method 
developed herein develops along analogous lines to the 
MCTDH approach and can be best thought of as an "or- 
bital" free approach since we work entirely at the level 



of the nuclear TV-body density. We do this by first writ- 
ing the configurational density n(l • • ■ N) that describes 
the statistical likelihood of finding the system in a given 
multi-dimensional configuration {r% ■ ■ ■ r^} as a superpo- 
sition of statistical approximates p{r\ ■ ■ ■ rjsr, c m ) that are 
joint probabilities for finding system at {r\ ■ ■ - rjv} and 
that it is a variant of some statistical distribution de- 
scribed by the approximate. These approximates can be 
any elementary probability distribution function that can 
be specified in terms of its statistical moments, c m , the 
simplest of which for our purposes are multi-dimensional 
gaussians. In this case, we need to be able to specify 
m(3iV(3iV + l)/2 + 3N + 1) = 0(mN 2 ) variables cor- 
responding to the elements of the covariance matrix, the 
central mean, and amplitude for m 3./V-dimensional gaus- 
sians. 

Explicit correlations between various degrees of free- 
dom can be excluded in straightforward way by factoring 
the approximates. For example, if we factorize the full co- 
variance matrix into individual atomic components, the 
configurational density can be cast in terms of the indi- 
vidual atomic densities 
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n(l ■ ■ ■ N) = Y[; 



(1) 



We can then expand each atomic density rii(i) as a linear 
combination of density approximates. 



ni(n) 
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(2) 



This dramatically reduces the number of coefficients we 
need to determine to mN x (6 + 3 + 1) = 0(mN). Inter- 
mediate factorization scheme yield similar scaling behav- 
ior allowing one to tune the computational complexity 
of the system depending upon the degree of correlation 
required by a particular physical problem. For example, 
one can define quasi-atoms by explicitly including covari- 
ance between the degrees of freedom of 2 or more atoms. 
As we shall derive next, each quasi-atom or atom will 
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then evolve in the mean-field of the other quasi-atoms of 
the system. 

In this paper, we present a grid-free adaptive hydro- 
dynamic approach for computing the quantum ground- 
state density for a system of N nuclei. Our approach uses 
Bayesian analysis to deduce from a statistical sampling 
of the density the best set of m statistical approximates 
describing that density. We then use a quantum hydro- 
dynamical scheme to move the sample points towards a 
minimal energy configuration. As proof of concept we 
consider the zero-point energy of small 4 and 5 atom 
clusters of Argon and Neon with pair-wise interatomic 
potential interactions. Finally, we discuss how the ap- 
proach may be used to develop new quantum-classical 
and fully quantum mechanical approaches for treating 
quantum mechanical solute particles (such as an excess 
e~ or He atom) in a liquid of classical or quasi-classical 
atoms (such as Ar or Ne). 



II. SELF-CONSISTENT FIELD EQUATIONS 

We begin by writing the full many-body Hamiltonian 
for the nuclear motion of a collection of atoms with pair- 
wise interaction potentials. 
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(3) 



where the first is the sum of the kinetic energies of the 
individual atoms and the second is the sum of the poten- 
tial energy contributions. For an arbitrary iV-body trial 
density, the energy functional is given by 

E[n}=T[n}+J2 [ / n^n^r^V^dndrj. (4) 

Since the kinetic energy operator is separable and we 
have assumed distinguishability amongst the constituent 
atoms, the kinetic energy term is the sum of the individ- 
ual kinetic energy functionals. 



As in electronic structure DFT, evaluating the kinetic 
energy functionals in an orbital free form is problematic 
since evaluating the quantum kinetic energy operator is a 
non-local operator and the density is a local function^, 

If instead we write the quantum wave function in polar 
form, as in the hydrodynamic formulation of quantum 
mechanic3^i£*i&, 



(6) 



T[n{l---N)] = Y,Ti[n{n)l 



(5) 



we can arrive at a stationary condition that if V0 = 0«£, 



V(l ■ ■ ■ N) - ^ J== V- V^R = const, (7) 

at all points in space. The constant is of course the en- 
ergy. By inspection, then, we can define our kinetic en- 
ergy functional as 



T[n(n)] = -2^7 / ^V^y^n)^. (8) 



Integrating by parts and letting n(i) — * at ±oo 
produces the familiar von Weizsacker kinetic energy 
functional^ 

1 f 1 - 

Tw[n(ri)] = +— / ——-\7in(n) • V '^{r^dr^ (9) 

oTTt J TL\Ti J 

Thus, the total energy functional is given in terms of the 
single particle densities. 



ni(ri)ni(rj)V(ij)dridrllO) 

i=i i±j J J 



Taking the variation of E[n] with respect to the 
single-particle densities with the constraint that 
Si/ ni{n)di = N, 



§ j f T w[ni(n)] + ^2 J J ni(ri)nj(rj)V(ij)dridrj - n yj m(ri)dri - 1 



= 0, (11) 



I 

leads to the following Euler-Lagrange equations: When satisfied, /x is the vibrational ground-state energy 

and the n^r^) = |</>i(i)| 2 are the probability densities of 
STwln^ri)} \ - f Tr ,..^ i \ j n /, x the individual nuclei. 

c / x + 2^ / V K l 3) n A r o) dr o - M = 0. (12) 

on n r t) y 4 J Let us take a simple pedagogic case of a particle in 
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a harmonic well, takin g the trial density density to be 
a Gaussian, n(x) — -\/a/7rexp(— ax 2 ). Evaluating the 
energy functional yields: 

1 ™,.,2 



E[n(x)\ = a - 

4m 



4a 



Minimizing with respect to the trial density 

SE[n] _dE _ Q 
Sn da 

yields the familiar E = cj/2 and a = mu>. This idea 
is easy to extend beyond purely harmonic systems and 
gaussian trial functions. Since n(i) is a probability dis- 
tribution function, it is a positive, real, and integrable 
function. 

In the next section, we show how the single particle 
densities can be estimated as super-positions of single 
particle density approximates and that the coefficients 
(rather moments) of the approximates can be optimized 
to compute both the ground state energy and the single- 
particle densities. 



III. MIXTURE MODELING 

The single-particle probability distribution functions 
(PDF) can be represented by a mixture modeli2i2£ by 
summing a finite number M of density approximates 



M 



n(r) = ^i>(r,i 



(13) 



where p(r, cm) is the probability that a randomly chosen 
member of the ensemble has the configuration r and is a 
variant of the mth approximate designated by c m . These 
approximates may be Gaussians or any other integrable 
multi-dimensional function which can be parameterized 
by its moments. For gaussian clusters, we have a weight 
p(c m ), a mean position vector /j, m , and a covariance ma- 
trix Cyyi • 

By definition, each joint probability in Eq. I13l is related 
to a pair of conditional probabilities according to the 
relation 

p(r, c m ) = p(c m )p{r\c m ) = n(r)p(c m \r), (14) 

where the forward conditional probability p(r\c m ) refers 
to the probability that a randomly chosen variant of c m 
has the configuration r. Conversely, the posterior prob- 
ability p(c m \r) refers to the probability that the config- 
uration point r is a variant of the approximate c m . In 
probability theory, n(r) and p(c m ) are known as marginal 
probabilities; however, we shall simply refer to them as 
the quantum density and weight of the mth approxi- 
mate, respectively. The expansion weights are strictly 
positive semidefinite and sum to unity. Substituting the 
first equality of Eq. [21 into Eq. ^5] we have 
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n(r) = ^2p(c m )p{r\c m ). 



(15) 



We have considerable freedom at this point in specify- 
ing the exact functional form of the conditional probabil- 
ities as well as the degree of correlation within each con- 
ditional. This freedom of specification allows us to con- 
struct "models" that explicitly take into account nonsep- 
arable correlations in configuration space. For the case 
of gaussian approximates this is accomplished by keeping 
or discarding various off-diagonal terms incorporates in 
the covariance matrix, C, 



P(r\c m ) 




(rd-/J m ,i).C m .(rd-At m ,d) 



(16) 



The term, || C -1 1|, stands for the reciprocal of the great- 
est value of the determinant of the covariance matrix. It 
is also possible to construct a model that assumes that 
each approximate is completely separable and takes the 
form of a product over the ^-dimensional configuration 
space, that reduces the covariance matrix, C, to a vari- 
ance vector, 



p(r|c TO ) = Y[ 



1 



■exp((r d - / x m , d ) 2 /(2a^ id )). (17) 



Numerical tests by Maddox and Bittner indicate that 
separable case is computationally faster for high dimen- 
sional systems, but produces a less accurate estimate of 
the quantum ground-state energy^i For larger systems 
the noncovariant case can certianly be used to speed 
calculations. We will examine the case where there is 
nonzero covariance between the three dimensions, but 
the atomic degrees of freedom have zero overlap. For a 
discussion of the strengths and weaknesses involved with 
mixture models, one is referred to the RefsiSLSi. This 
approximation provides a sufficient approximation to the 
density for the calculation of the ground state. 

Once a model is decided upon one must then deter- 
mine the parameters, in this case the Gaussian parame- 
ters p(c m ), p mi and C' m , for each approximate from the 
statistical points representing the density. For instance 
the mean position vectors of the approximates are defined 
by the moments of the forward conditional probabilities 



Mm = J rp(r\c m )dr. 



(18) 



Rearranging Eq. ^]and substituting into Eq. ^[ we can 
write these parameters as 



n(r)p(c m r) 

Pm = / r — tir, 

p{c m ) 



(19) 



this is easily approximated by summing over an ensemble 
of points {r„} sampled from the n(r) PDF, 



^ m ~ l^TT — \ r n p{c m \r n ). 



(20) 
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We also define similar expressions for the covariance ma- 
trix and the expansion weights 



1 N 



(21) 



Taking the variation of L with respect to the model pa- 
rameters generates a series of update-rules for moving the 
approximates through parameter space in the direction 
along V Cm Lm&. For the case of Gaussian approximates, 
the update rules for the mean, covariance matrix, and 
marginal probabilities are given by, 



1 



N 



Np(c m ) 



^(r„ - fi) T (r„ - n)p{Cn 



(22) 



For the separable case, the variances are given by the 
diagonal elements cr^ i = (C m )u. The posterior terms 
p(c m \r n ) for each r n sample point in Eqs. 1201221 are eval- 
uated directly from the forward probabilities according 
to Bayes' equation, 



p{c m \r n ) 



p(c m )p(r n \c m ) 

EmP( c m)p( r ™l c m)' 



(23) 



Within this viewpoint, the sample points can be con- 
sidered to be a data set that represents the results of a se- 
ries of successive measurements. Each data point carries 
an equal amount of information describing the underlying 
quantum probability distribution function. Bayes' equa- 
tion gives the ratio of how well a given estimate describes 
r„ to how well r n is described by all of the approximates. 
Thus, it represents the fraction of explanatory informa- 
tion that a given sample point gives to the m-th approx- 
imate. The estimate which best describes the point will 
have the largest posterior probability at that point. Eqs. 
I2UI23I can be iterated self-consistently in order to deter- 
mine the best possible set of parameters that describe 
n(r) in terms of a given ensemble of data points. In 
doing so, we effectively maximize the log-likelihood that 
the overall density model describes the entire collection 
of data points. 



L = log]Jn(r„) 



~ m V7 t 



<5C„, — 



2(C m C m ) , 
Np{c m ) 



. L, 



Sp{c m ) = — (diagffi] 



Where ® is the Kronecker product, is the vector of all 
expansion weights, = [p(ci), . . . ,p(c m )] T , and diag[0] 
is a matrix with the elements from f2 constituting the 
diagonal entries^ 

The expectation maximization algorithm described 
above allows us to generate an approximate analytical 
functional form for the single particle density via sta- 
tistical sampling over an ensemble of points. The next 
step is to adjust the single-particle densities themselves 
to produce a lower total energy. We do this by de- 
riving the quantum hydrodynamic equations of motion 
for the sample points. The quantum Hamilton- Jacobi 
equation generates the equations of motion for the ray- 
lines of a time-dependent solution to the Schrodinger 
equation, 24 ! 25 ! 26 ! 27 . 



dS x 



2m,- 



Yl ff niirijUjir^Viifidridrj - ^ 



iV.V^fc) = 



(24) 



r 



Since the density is separable into components, we eas- 
ily arrive at a set of time-dependent self-consistent field 
equations whereby the motion of atom i is determined 
by the average potential interaction between atom i and 
the rest of the atoms in the system. 



V?V^) = 0- (25) 



2rrii 
1 



2m 



Taking ViS = p = m^r as a momentum of a particle, 
the equations of motion along a given ray-line or char- 
acteristic r n (t) of the quantum wave function are given 
by 

miVn = (^(y)) n i( r i)*j - ViQ[ni(ri)] (26) 

where Q[n(i)] is the Bohmian quantum potential speci- 
fied by the last term in Eq. [23 Stationary solutions of 
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time-dependent Schrodinger equation are obtained when- 
ever m,r n = 0. Consequently, by relaxing the sample 
points in a direction along the energy gradient specified 
by 




keeping n(rj) fixed. This generates a new statistical sam- 
pling, which we then use to determine a new set approx- 
imates. 

This process is similar to the semiclassical approxima- 
tion strategy for including quantum effects into other- 
wise classical calculations introduced by Garaschuk and 
Rassolo\*22i2S. This semiclassical approximate method- 
ology is based upon de Broglie-Bohm trajectories and 
involves the convolution of the quantum density with 
a minimum uncertainty wave packet which is then ex- 
panded in a linear combination of Gaussian functions 

p(x) « f(x) =Y, c l exp[-<4(z - X n )}. (28) 
n 

The Gaussian parameters s — {c„, X n , a n } in Eq. 1281 are 
determined by minimizing the functional 

F = J [p(x) - f(x)] 2 dx (29) 

using an iterative procedure which explicitly involves 
solving the set of equations dF/dsk = 0. The 
parametrized density leads to an approximate quantum 
potential (AQP) that is used to propagate an ensemble 
of trajectories. This approach has been used successfully 
in computing reactive scattering cross-sections for the co- 
linear H+H 2 reaction in one dimension. 

It is important at this point to recognize the numer- 
ical difficulties our group and others have faced in de- 
veloping hydrodynamic trajectory based approaches for 
time-dependent systemsi 21 i 27 i 30 i 31 i 32 i 33 The foremost dif- 
ficulty is in the accurate evaluation of the quantum po- 
tential from an irregular mesh of points^S* 3 ^ The quan- 
tum potential is a function of the local curvature of the 
density and can become singular and rapidly varying as 
nodes form in the wave function or when the wave func- 
tion is sharply peaked, i.e. when n 1 / 2 — > faster than 
V 2 ?! 1 / 2 — * 0. These inherent properties make an accu- 
rate numerical calculation of the quantum potential and 
its derivatives very difficult for all but the simplest sys- 
tems. These difficulties are avioded in the cluster model 
approximation of the density, using the expectation max- 
imization (EM) algorithm^. By obtaining a global op- 
timal function that describes the density, we can analyt- 
ically compute the quantum force with great accuracy. 
The issue of nodes is essentially avoided so long as we 
are judicious in our choice of density approximates. If we 
choose node-free approximates, then our overall density 
will likewise be node free. For the purpose of determining 
vibrational ground-states, this seems to be a worthwhile 
compromise. 

The algorithm can be summarized as follows 



1. For each atom, generate and sample a normalized 
trial density riifri). 

2. Using the EM routines and the given sample of 
points, compute the coefficients for the density ap- 
proximates. 

3. Compute the forces on each point using Eg. 1771 and 
advance each point along the energy gradient for 
one "time" step, either discarding or dampening 
the velocity of each point. This generates a new 
sample of points describing the single-particle den- 
sity for each atom. The new distribution should 
have a lower total energy since we moved the sam- 
ple points in the direction towards lower energy. 

Iterating through these last two steps, we rapidly con- 
verge towards the energy minimum of the system. 



IV. VIBRATIONAL GROUND STATE OF 
RARE-GAS CLUSTERS 

As proof of concept we examine the ground vibrational 
state energies of Lennard- Jones clusters. The simple 
Lennard-Jones pair potential provides reliable informa- 
tion about complicated systems, and has been used in a 
number of recent studies with just a few examples listed 
in Refsi 3 ^ 3 ^* 3 ^. Ground state energies of rare gas clus- 
ters are easily enough modeled with molecular dynamics 
simulations, but the quantum corrections are often quite 
large for these cases. These corrections are important be- 
cause the quantum character strongly affects the thermo- 
dynamics via changes in the ground state structure due 
to increasing zero-point energies^. The zero point energy 
corrections for the small clusters modeled here can be up 
to 0.66kJ/mole. Indeed, quantum corrections have been 
shown to lower solid to liquid transition temperatures by 
approximately 10%, and the zero point energy for small 
clusters can account for up to 35% of the classical binding 
energy^. 

The effects from quantum derealization are intuitively 
understood in the present approach through the quantum 
potential term in the equations of motion. This explains 
why the quantum derealization can account for a low- 
ering of the transition temperature because some kinetic 
energy is always present even at very low temperatures. 
This spreading of the wave packet is known as a "soft- 
ening" of the crystal which leads to a lowering of the 
melting temperaturei 3 ^ These effects have been studied 
in the context of the transition from molecular to bulk- 
like properties of clusters. 

In the calculations presented here, we used 300 sta- 
tistical points to represent the density of each atom 
and propagated the SCF equations described above until 
the energy and the density were sufficiently converged. 
Typically, this required 1.5 million to 3 million cycles. 
Along the course of the energy minimization, we strongly 
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FIG. 2: The average potential energy (V) and total energy 
(Q) + (V) of the Ars and Nes clusters in kj/mole. The steps 
are measured in millions. 




FIG. 1: The isodensity contour plots of the clusters at a value 
of 0.006. In the upper left is the Ar4 cluster, in the upper right 
is the Ne4, lower left has the Ars, and then bottom right is 
Nes. The axis are listed in atomic units. 



TABLE I: Inter-atomic distances for X5 clusters in atomic 
units. 



distances 


Argon 


Argon (cl) 


Neon 


Neon (cl) 


rdi,2 


7.343±.019 


7.225 


6.167±.045 


5.927 


rdi,3 


7.327±.018 


7.225 


6.115±.031 


5.927 


rdi,4 


7.288±.018 


7.199 


6.057±.035 


5.906 


rdi,5 


7.269±.016 


7.199 


6.048±.027 


5.906 


rd-2,3 


7.339±.023 


7.225 


6.162±.055 


5.927 


rd 2 ,A 


7.276±.018 


7.199 


6.060±.039 


5.906 


rd 2 , 5 


7.285±.016 


7.199 


6.070±.035 


5.906 


rd-j.,4, 


7.266±.018 


7.199 


6.112±.060 


5.906 


rd 3 , 5 


7.271±.016 


7.199 


6.077±.035 


5.906 


rd 4 , 5 


11.832±.020 


11.7349 


9.848±.039 


9.627 



damped the time-evolution of the sample points to elim- 
inate as much of the oscillations and breathing of the 
density components as possible. 

The Lennard- Jones parameters for the Argon clusters 
are e = 0.9976 kJ/mole and a = 3.421, and e = 0.3059 
kJ/mole and a = 2.79 A for the Neon clusters^. Initial 
configurations for the simulations are chosen to be close 
to the classical molecular dynamics minimum energy ge- 
ometry, and are given some initial Gaussian spread. 

We show in Fig. ^ isodensity (0.006) contour plots 
for the Ai'4, Ars, Ne4, and Nes4i. One can see quite 
clearly the underlying three-dimensional shape of the 
cluster along with the derealization of each atom about 
its central location. Each density "lobe" is nearly spher- 
ical with some elongation. These density plots give a 
suggestive view of the overlap of the densities which is 
ignored in Eq. ^and therefor ultimately in Eq. For 
this system this overlap turns out to be minor, but for 
atoms such as Helium this would have to be taken into 
account. 

In Fig. we show the total energy an the total poten- 
tial energy for the Ars and Nes clusters as the system 
converges towards its lowest energy state. Initially there 
is a rapid restructuring of the densities as they adjust 
to find a close approximation to the actual ground state 
density. Following this initial rapid convergence, there 
is slower convergence phase as the density is further re- 
fined. During this process as the sample points look for 
a configuration which fully equalizes of the quantum and 
kinetic energy terms from Eq. 0the density approxima- 



tion can sometimes prove inadequate and points can be 
pushed into temporary higher energy regions. This leads 
to the fluctuations seen in the energies and any other 
averaged quantity, such as the interatomic distances. In 
order to compute meaningful values for the energy and 
distances, we averaged these quantities over the last half 
million or so cycles. As can be seen from Figure [21 the 
Nes cluster is slower to converge, but eventually does so 
around step 200 (million). 

Tables U and HJ lists the averaged interatom distances 
for each cluster compared to the equilibrium distances for 
the corresponding classical case. For the case of Ars the 
numerical fluctuations lead to an uncertainly of about 
0.3% in the interatomic distances and for -/Ves, a 0.5% 
uncertainty in the interatomic distances, it is important 
to note that the fluctuations mentioned here get smaller 
for larger systems as can seen by comparing the results 
for Ars with Ar4. This is because the well depths for the 
sample points become more pronounced. This has impor- 
tant implications since we hope to extent this method to 
larger systems. Ne4 can be seen to have the largest fluc- 
tuations, which is expected since it is the most quantum 
mechanical. All in all these values compare well with the 
classical distances. In general, the quantum distances 
are slightly larger due to the fact that the gaussian atom 
densities are sampling part of the anharmonic attractive 
portion of the pair-potential. 

Table Hill summarizes the various contributions to the 
the total energy for each cluster . The "classical" ener- 
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TABLE II: Inter-atomic distances for X4 clusters in atomic 
units. 



distances 


Argon 


Argon fell 


Neon 


Neon (cl) 


rd\,2 


7.362±.024 


7.21 


6.167± .043 


5.918 


rdi,3 


7.296±.017 


7.21 


6.104 ± .051 


5.918 


rd 1A 


7.294±.019 


7.21 


6.092 ± .049 


5.918 


rd2,3 


7.320±.022 


7.21 


6.140 ± .038 


5.918 


rd 2 ,4 


7.305±.016 


7.21 


6.077 ± .027 


5.918 


rd 3A 


7.286±.017 


7.21 


6.094 ± .040 


5.918 



gies, V c are the energy minimum of the total potential 
energy surface corresponding to the classical equilibrium 
configuration. (V), (Q), and (E) are the total quantum 
potential energy, kinetic energy, and total energy of each 
cluster. The difference between the classical potential 
minimum, V c and (E) is the zero point energy. For our 
results we calculate a virial-like term which is the ra- 
tio of the kinetic to the potential energy of the system, 
(Q)/(V). It turns out the values we get for this term, 
seen in Table ITTT1 are remarkably close to the values of 
the de Boer parameter, A, for the atoms we are examin- 
ing. 



TABLE III: Table of the energies of the clusters in kj/mole. 





A1-4 


Ne 4 


Ar 5 


Ne 5 


Vc 


-11.972 


-3.655 


-18.166 


-5.545 


(V) 


-11.337±.225 


-3.184± .056 


-17.302 ± .344 


-4.895±.055 


(Q) 


0.462±.033 


0.460 ± .034 


0.669 ± .041 


0.646±.036 


(E) 


-10.874±.215 


-2.724 ± .049 


-16.632 ± .330 


-4.249±.043 


(Q)/(V) 


0.041 


0.144 


0.038 


0.132 



The de Boer parameter has been used in attempts 
to rationalize the effects of quantum derealization for 
Lennard-Jones systems AiSi^i Each Lennard- Jones sys- 
tem can be defined in terms of its parameters e, or its 
potential depth, its length scale of a, and mass m. For 
a given set of parameters, the thermal de Broglie wave- 
length, A = ft 1 'a 'v 'niksT provides a means of approx- 
imating the quantum effects, at some reduced temper- 
ature of the system, T* — ksT/e. Further taking the 
ratio of A for two different sets of parameters provides a 
means of comparing the quantum effects of one system 
versus the other. This leads to the de Boer parame- 
ter, A = h/r m y/me, which is the ratio of the de Broglie 
wavelength of an atom with energy e, with an intermolec- 
ular distance, r rn . Basically, the de Boer parameter is 
useful for comparing the quantum character of different 
Lennard-Jones clusters or liquids at a given temperature, 
in our case zero4& A has a classical limit of A = and 
anything above A 0.3 is considered a quantum system. 
For Argon the de Boer parameter is, A 0.03 which cor- 
responds to a classical system, and quantum effects can 
be treated as a perturbation. For Neon the de Boer pa- 
rameter turns out to be A « 0.1 which is a quasi- quantum 
system. 

The de Boer parameter measures the derealization of 
the system compared to its size. The virial like term 
measures the percentage energy contained in the kinetic 
term of the Hamiltonian. The kinetic energy, also called 
the quantum potential energy, is also a measure of the 



derealization of the system. So we see that the two terms 
are essentially measuring the same entity. For our results 
we also see a possible trend to smaller values of (Q)/(V), 
as the system gets larger. 



V. DISCUSSION 

A method for calculating ground state configuration 
of quantum clusters and liquids has been outlined based 
upon some previous work in approximating densities as 
quantum statistical distribution function. The quantum 
and the Lennard-Jones potentials are used to propagate 
an ensemble of Monte Carlo statistical points, in a DFT 
like procedure. This is an orbital free approach in the 
sense that we only work at the level of the nuclear den- 
sity. In order to do this we outline a "cluster" model and 
expectation maximization (EM) algorithm which is used 
to obtain the density in terms of the statistical points 
representing each atom. The Lennard-Jones potential is 
calculated in a mean field sense by averaging over the sta- 
tistical points of each atom, and the quantum potential 
is calculated from the density obtained in the EM algo- 
rithm. Results were presented for 4 and 5 atom clusters 
of Argon and Neon. The results indicate good agreement 
with the general classical results, but the quantum cor- 
rections can be seen to be significant. Also shown is that 
the virial term we measure to approximate the quantum 
effects is related to the de Boer parameter used in previ- 
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ous studies. 

The method outlined also seems to provide a means of 
artificial control of the amount of quantum mechanical 
information desired from a calculation. The covariances 
between atoms may be set to zero as we have done, which 
corresponds to distinguishable particles, but keeping the 
interatomic covariances seems to provide a possible path 
to including other effects such as exchange energies and 
the like. This all comes at the expense of increased com- 
plexity and reduced computational speed. 

One could also tune the amount of quantum effects by 
making ft a parameter that can vary between the values 
of to 1, in atomic units, in the equation for the quantum 
potential. Additionally the method is easily extended to 



potentials beyond the Lennard-Jones. For instance, co- 
valent bonds could be modeled with harmonic oscillators 
or Morse potentials, and Coulomb potentials are easily 
modeled. 
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